สถิติแบบ Frequentist

สถิติแบบดั้งเดิมหรือที่เรียกกันว่า frequentist เป็นสถิติวิเคราะห์ที่ถูกใช้งานเป็นหลักในปัจจุบัน สถิติเกือบทุกตัวที่ผู้เรียนได้เรียนมาจนถึงปัจจุบันนี้ล้วนเป็นสถิติวิเคราะห์ภายใต้กระบวนทัศน์แบบ frequentist ทั้งสิ้น สถิติแบบ frequentist นี้มีข้อตกลงเบื้องต้นว่า (1) พารามิเตอร์เป็นค่าคงที่ ที่มีอยู่จริงในประชากรเป้าหมาย อย่างไรก็ตามด้วยข้อจำกัดที่ไม่สามารถเก็บรวบรวมข้อมูลจากทั้งประชากรได้ทำให้ผู้วิเคราะห์ไม่ทราบค่าพารามิเตอร์ดังกล่าว และจะประมาณโดยใช้ข้อมูลจากตัวอย่าง (2) ข้อมูลตัวอย่าง (sample data) เขียนแทนด้วย \(\bf{x}= \left \{\bf{x_1}, \bf{x_2},...,\bf{x_n} \right \}\) เป็นค่าที่สุ่มมาจากการประชากร ดังนั้นตัวอย่างที่สุ่มได้จึงเป็นค่าสุ่มที่และมีความเป็นไปได้ทั้งหมดเท่ากับ เมื่อ \(N\) คือขนาดประชากร และ \(n\) คือขนาดตัวอย่าง จากข้อตกลงเบื้องต้นข้อที่ (2) นี้จึงทำให้ (3) ค่าสถิติที่ประมาณได้จากข้อมูลตัวอย่างเป็นตัวแปรสุ่มที่มีการแจกแจงความน่าจะเป็น เรียกการแจกแจงความน่าจะเป็นของค่าสถิตินี้ว่า การแจกแจงความน่าจะเป็นของฟังก์ชันที่ได้จากตัวอย่างสุ่ม (sampling distribution)

ลักษณะการแจกแจงความน่าจะเป็นของค่าสถิติดังกล่าวมักประมาณได้ด้วย 2 วิธีการ วิธีการแรก ประมาณโดยใช้ทฤษฎีลิมิตลู่เข้าสู่ส่วนกลาง (central limit theorem) ซึ่งเป็นการคาดการณ์ลักษณะการแจกแจงของค่าสถิติในเชิงทฤษฎี เช่น การประมาณช่วงความเชื่อมั่นของค่าเฉลี่ยมีสูตรเป็น \(\overline{X}\pm t_{\alpha/2,n-1} \frac{S}{\sqrt{n}}\) จะเห็นว่าการประมาณช่วงความเชื่อมั่นดังกล่าวมีการใช้การแจกแจงที (student t distribution) ในการคำนวณส่วน margin of error การแจกแจงดังกล่าวพิสูจน์ในทางทฤษฎีโดยใช้ central limit theorem ได้ว่า เป็น sampling distribution ของ \(\overline{X}\) หากข้อมูลตัวอย่างที่ใช้ในการวิเคราะห์สุ่มมาจากประชากรที่มีการแจกแจงแบบปกติหรือใกล้เคียง แต่ไม่ทราบค่าความแปรปรวนของประชากร หรือถ้าหากประชากรไม่ได้มีการแจกแจงแบบปกติ แต่ตัวอย่างมีขนาดใหญ่เพียงพอ ก็ยังสามารถใช้ sampling distribution ดังกล่าวประมาณการแจกแจงของ \(\overline{X}\) ได้อยู่ หรือในการทดสอบสมมุติฐาน \(H_0: \mu=\mu_0\) สถิติทดสอบที่ใช้คือ \(t^*=\frac{\overline{X}-\mu_0}{S/ \sqrt{n}}\) จะมี sampling distributionในทางทฤษฎีเป็นการแจกแจงแบบที ที่มีองศาความเป็นอิสระเท่ากับ \(n-1\) เช่นเดียวกันหากข้อกำหนดเบื้องต้นของ central limit theorem เป็นจริงหรือใกล้เคียง การทราบการแจกแจงของค่าสถิติดังกล่าวทำให้ผู้วิเคราะห์สามารถประเมินความคลาดเคลื่อนในการตัดสินใจของการทดสอบสมมุติฐาน และสรุปผลได้ดังที่ได้เรียนไปแล้วในรายวิชาพื้นฐาน

ความถูกต้องของทฤษฎีนี้ขึ้นอยู่กับปัจจัยแวดล้อมที่เกี่ยวข้องว่ามีความสอดคล้องกับทฤษฎีมากน้อยเพียงใด ปัจจัยหนึ่งที่มีความสำคัญคือขนาดตัวอย่าง กล่าวคือหากขนาดตัวอย่างมีค่าน้อยเกินไป หรือการแจกแจงของข้อมูลมีความเบี่ยงเบนออกไปจากข้อตกลงเบื้องต้นของการวิเคราะห์ไปมาก ค่าสถิติที่คำนวณได้อาจไม่มีลักษณะการแจกแจงที่เป็นไปตามทฤษฎี และส่งผลให้การอนุมานเชิงสถิติมีความคลาดเคลื่อน อีกวิธีการหนึ่งที่สามารถใช้ในการประมาณ sampling distribution ของค่าสถิติได้คือใช้วิธีการสุ่มซ้ำ (resampling method) เช่น boostraping หรือ jackknife เพื่อประมาณแนวโน้มการแจกแจงของค่าสถิติ เป็นต้น

t.test(dat$Ach, mu=5.5)

## 
##  One Sample t-test
## 
## data:  dat$Ach
## t = 5.078, df = 149, p-value = 1.123e-06
## alternative hypothesis: true mean is not equal to 5.5
## 95 percent confidence interval:
##  5.709732 5.976934
## sample estimates:
## mean of x 
##  5.843333

คำถาม

1. output ข้างต้นแสดงผลการประมาณช่วงความเชื่อมั่น และทดสอบสมมุติฐานทางสถิติเกี่ยวกับผลสัมฤทธิ์ทางการเรียนของนักเรียน (คะแนนเต็ม 10) โดยจากการประมาณด้วยช่วงความเชื่อมั่น 95% ในข้างต้นพบว่ามีค่าเท่ากับ \([5.71, 5.98]\) ซึ่งหมายความว่าอะไร?

2. จากผลการทดสอบสมมุติฐาน \(H_0: \mu_{Ach}=5.5\) พบว่ามีค่า p-value < .000 หมายความว่าอะไร?

3. เราสามารถตัดสินใจยอมรับสมมุติฐาน \(H_0\) ได้หรือไม่?

จากข้อตกลงเบื้องต้นของสถิติแบบ frequentist นี้ทำให้การวิเคราะห์ข้อมูลเกิดข้อจำกัดต่าง ๆ ดังนี้

1. ด้านที่สำคัญคือด้านการแปลผลการวิเคราะห์ที่จะเห็นว่าความหมายของผลการวิเคราะห์ต่าง ๆ เป็นการกล่าวถึงพารามิเตอร์ที่สนใจแบบอ้อม ๆ ทั้งสิ้น ไม่มีการวิเคราะห์ใดที่สามารถอ้างอิงไปยังพารามิเตอร์ที่สนใจได้โดยตรง ทั้งนี้เป็นเพราะการอนุมานเชิงสถิติแบบ frequentist พึ่งพาเครื่องมือหลักคือ samping distribution ของค่าสถิติ ในขณะที่พารามิเตอร์ถูกสมมุติให้เป็นค่าคงที่ตั้งแต่ต้น probability statement ต่าง ๆ จึงเป็นของค่าสถิติทั้งสิ้น ดังความความเชื่อมั่นที่รับประกันในช่วงความเชือมั่น และค่าความคลาดเคลื่อนในการทดสอบสมมุติฐานจึงเป็นค่าความน่าจะเป็นที่ใช้รับประกันความเป็นไปได้ของค่าสถิติทั้งสิ้นไม่ใช่ของค่าพารามิเตอร์ที่สนใจจริง ๆ ในการวิเคราะห์

2. ด้านที่สองคือด้านการประมาณค่าพารามิเตอร์ภายใต้โมเดลซับซ้อน เครื่องมือสำคัญสำหรับประมาณค่าพารามิเตอร์ในโมเดลต่าง ๆ ของสถิติแบบ frequentist คือตัวประมาณ เช่น ตัวประมาณกำลังสองน้อยสุด ตัวประมาณภาวะความควรจะเป็นสูงสุด ซึ่งล้วนเป็นวิธีการที่จะหาค่าประมาณที่ดีที่สุดโดยอิงจากจุดต่ำสุด หรือจุดสูงสุดของฟังก์ชันวัตถุประสงค์ ดังตัวอย่างในรูปด้านล่าง ซึ่งในกรณีที่โมเดลมีความซับซ้อนมาก ๆ ฟังก์ชันวัตถุประสงค์ดังกล่าวจะมีความซับซ้อนโดยมีจำนวนมิติที่สูงขึ้น และอาจมีจุดสูงต่ำ และ/หรือ จุดต่ำสุดสัมพัทธ์จำนวนมาก ซึ่งเป็นอุปสรรคให้การประมาณค่าด้วยวิธีการในข้างต้นอาจได้ค่าประมาณที่คลาดเคลื่อนไปจากค่าที่เหมาะสมที่สุด หรืออัลกอริทึมไม่สามารถหาค่าประมาณที่เหมาะสมได้ภายใต้จำนวนรอบของการประมาณที่กำหนด

3. สืบเนื่องจากที่การอนุมานเชิงสถิติแบบ frequentist นั้นมักอิงกับทฤษฎีความน่าจะเป็นที่อิงข้อกำหนดเบื้องต้นที่ขนาดตัวอย่างต้องมีขนาดใหญ่เพียงพอ ซึ่งทำให้สถิติแบบ frequentist มักทำงานได้ไม่ดีในกรณีที่ตัวอย่างมีขนาดเล็ก

สถิติแบบเบส์ (Bayesian Statistics)

สถิติแบบเบส์เป็นกระบวนทัศน์ในการวิเคราะห์ข้อมูลอีกลักษณะหนึ่งที่มีความแตกต่างจากไปกระบวนทัศน์แบบดั้งเดิมที่นักวิเคราะห์ข้อมูลทั่วไปคุ้นเคยกัน ซึ่งทำให้สถิติแบบเบส์มีจุดเด่นและสามารถแก้ไขข้อจำกัดที่พบในสถิติแบบ frequentist ได้ ความแตกต่างดังกล่าวเริ่มตั้งแต่ข้อตกลงเบื้องต้นของสถิติแบบเบส์ที่แตกต่างจากสถิติแบบดั้งเดิมอย่างมาก กล่าวคือ

(1) พารามิเตอร์เป็นตัวแปรสุ่ม สถิติแบบเบส์มองว่าพารามิเตอร์เป็นค่าที่ผู้วิเคราะห์สนใจจะศึกษาและต้องการทราบอย่างแท้จริง แต่อย่างไรก็ตามเนื่องจากผู้วิเคราะห์ไม่ทราบข้อมูลประชากร การคาดการณ์หรือความเชื่อ (belief) เกี่ยวกับพารามิเตอร์ดังกล่าวจึงมีความไม่แน่นอน สภาพดังกล่าวจึงสมเหตุสมผลที่จะกำหนด probability statement ให้กับพารามิเตอร์ในโมเดล

จากข้อกำหนดนี้ผลการประมาณค่าพารามิเตอร์ที่ได้จากวิธีการแบบเบส์จึงจะไม่ได้ค่าประมาณแบบค่าเดียว แต่จะได้เป็นการประมาณการแจกแจงของพารามิเตอร์แทน เช่น การประมาณค่าเฉลี่ยผลสัมฤทธิ์ของนักเรียน ผลการวิเคราะห์แบบเบส์จะไม่ได้ให้เป็นค่าประมาณเพียงค่าเดียวคือ 5.84 คะแนน ดัง output ข้างต้น แต่จะให้เป็นการแจกแจงของพารามิเตอร์ค่าเฉลี่ยผลสัมฤทธิ์ทางการเรียน ซึ่งเรียกว่าการแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) และสารสนเทศเกี่ยวกับค่าเฉลี่ยของผลสัมฤทธิ์ดังกล่าวจะถูกสกัดออกมาจากการแจกแจงนี้

(2) ข้อมูลตัวอย่างเป็นค่าคงที่ ถึงแม้ว่าในความเป็นจริงตัวอย่างที่สุ่มมาจากประชากรจะมีความเป็นไปได้ที่หลากหลาย และมีความไม่แน่นอนในการได้มาซึ่งตัวอย่างแต่ละชุด อย่างไรก็ตามเมื่อดำเนินการสุ่มตัวอย่างมาทำการวิเคราะห์แล้ว ความเป็นไปได้จำนวนมากก่อนการสุ่มตัวอย่างจะเหลือเพียงความเป็นไปได้เดียว สถิติแบบเบส์จึงถือว่าข้อมูลตัวอย่างเป็นค่าคงที่ ไม่มีการกำหนด probability statetment ให้กับข้อมูลตัวอย่างนี้

สถิติแบบเบส์สามารถนำมาใช้เพื่อแก้ไขข้อจำกัดที่เกิดขึ้นจากสถิติแบบ frequentist ได้ดังนี้

1. จากข้อกำหนดเบื้องต้นข้างต้น การอนุมานเชิงสถิติด้วยสถิติแบบเบส์จะดำเนินการผ่านการแจกแจงความน่าจะเป็นของพารามิเตอร์โดยตรงที่เรียกว่า การแจกแจงความน่าจะเป็นภายหลัง (posterior distribution) การแจกแจงดังกล่าวเป็นเครื่องมือสำคัญในการอนุมานเชิงสถิติแบบเบส์ แทน sampling distribution ในการอนุมานเชิงสถิติแบบดั้งเดิม ซึ่งจุดนี้ทำให้สถิติแบบเบส์สร้างความแตกต่างอย่างมากเมื่อเปรียบเทียบกันสถิติแบบดั้งเดิม โดยจะเห็นว่าการประมาณค่าพารามิเตอร์แบบเบส์ไม่ได้อิงกับค่าสถิติเพียงค่าเดียวเหมือนกับสถิติแบบดั้งเดิม แต่อิงกับการแจกแจงความน่าจะเป็นภายหลัง

2. การที่สถิติแบบเบส์ใช้การแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์เป็นเครื่องมือหลักในการอนุมานเชิงสถิติ จึงทำให้การวิเคราะห์ทำได้อย่างตรงไปตรงมาโดยไม่ต้องอาศัยทฤษฎีหรือวิธีการที่ซับซ้อนเหมือนกับสถิติแบบดั้งเดิม กล่าวคือผู้วิเคราะห์ใช้เพียงสถิติพื้นฐาน เช่น ค่าเฉลี่ย มัธยฐาน ฐานนิยม ส่วนเบี่ยงเบนมาตรฐาน หรือเปอร์เซ็นไทล์ เป็นเครื่องมือในการสรุปสารสนเทศเกี่ยวกับพารามิเตอร์ที่สนใจออกมาจากการแจกแจงความน่าจะเป็นภายหลังที่ประมาณได้ ก็เพียงพอแล้ว นอกจากนี้การแปลผลเพื่ออนุมานเกี่ยวกับพารามิเตอร์ยังทำได้ง่ายและตรงไปตรงมา และได้ข้อสรุปโดยตรงไปยังปริภูมิของพารามิเตอร์ที่เป็นเป้าหมาย แตกต่างจากข้อสรุปที่ได้จากสถิติแบบดั้งเดิมที่เป็นข้อสรุปเกี่ยวกับการแจกแจงของค่าสถิติซึ่งอยู่บนปริภูมิตัวอย่างเท่านั้น

3. การประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ในโมเดล จะใช้สารสนเทศสองส่วนด้วยกัน ส่วนแรกเป็นสมมุติฐานหรือความเชื่อเบื้องต้น (prior belief) เกี่ยวกับพารามิเตอร์ที่สนใจ ซึ่งกำหนดอยู่ในรูปการแจกแจงความน่าจะเป็นที่เรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์ (prior distribution) และส่วนที่สองคือสารสนเทศที่ได้จากข้อมูลตัวอย่าง ซึ่งแตกต่างจากสถิติแบบดั้งเดิมที่จะใช้สารสนเทศจากข้อมูลตัวอย่างเท่านั้น

ลองพิจารณาตัวอย่างต่อไปนี้

สมมุติว่าต้องการคาดการณ์คะแนนสอบของนายบูม ที่มีความเป็นไปได้ 4 แบบดังรูป

ก่อนการเก็บรวบรวมข้อมูลจริงในชั้นเรียนของนายบูม สมมุติว่าผู้วิเคราะห์มีข้อมูลในอดีตที่เกี่ยวข้องกับคะแนนสอบของนายบูมดังนี้

1. ไม่มีข้อมูลใด ๆ ที่เป็นประโยชน์ (ให้น้ำหนักกับความเป็นไปได้ A1, A2, A3 และ A4 อย่างไร?)

2. หากทราบเพิ่มเติมว่าการสอบที่ผ่าน ๆ มาโดยเฉลี่ยนายบูมสอบได้คะแนนคิดเป็นอันดับที่ 4 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 จะเหมือนเดิมมั้ย?)

3. หากทราบเพิ่มเติมอีกว่า การสอบที่ผ่าน ๆ มา แก้วซึ่งเป็นเพื่อนในห้องเรียนของนายบูม สอบได้คะแนน 38 คะแนน คิดเป็นคะแนนที่น้อยที่สุดในชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

4. หากทราบอีกว่า ในการสอบที่ผ่าน ๆ มา นิดสอบได้คะแนนโดยเฉลี่ยคิดเป็น 89 คะแนน ซึ่งเป็นอันดับที่ 3 ของชั้นเรียนนี้ (น้ำหนักของความเป็นไปได้ทั้ง 4 เหมือนเดิมหรือไม่?)

การแจกแจงก่อนหน้าของคะแนนสอบนายบูม

ผู้วิเคราะห์สามารถปรับสารสนเทศในการแจกแจงความน่าจะเป็นก่อนหน้าในข้างต้นได้ หากมีการเก็บรวบรวมข้อมูลตัวอย่างเพิ่มเติม เช่น

สมมุติว่าการแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมเป็นแบบ uniform และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

อีกตัวอย่างหนึ่ง สมมุติว่า การแจกแจงความน่าจะเป็นก่อนหน้าของคะแนนนายบูมมีน้ำหนักสูงที่ A3 และมีข้อมูลการสอบครั้งล่าสุดเพิ่มเติมดังรูป

กำหนดตัวแปรในการวิเคราะห์

สมมุติว่าผู้วิจัยต้องการวิเคราะห์ความลำเอียงในการออกหน้าหัวหรือก้อยของเหรียญที่ผลิตขึ้นมารุ่นหนึ่ง การวิเคราะห์นี้ตัวแปรที่สนใจจึงเป็นผลลัพธ์ที่สังเกตได้จากการโยนเหรียญในแต่ละครั้ง กำหนดให้ \(y_i\) แทนค่าสังเกตที่ได้จากการโยนเหรียญในครั้งที่ \(i\) จะได้ว่า

\(y_i=\left\{\begin{matrix} 1 & ; H\\ 0 & ; T \end{matrix}\right.\)

กำหนดโมเดลของค่าสังเกต

เมื่อกำหนดตัวแปรที่สนใจแล้วจะพบว่าตัวแปรที่สนใจมีเพียงตัวแปรตาม 1 ตัว โดยมีค่าสังเกตที่เป็นไปได้เพียง 2 ค่าได้แก่ H หรือ T ซึ่งแทนด้วย 1 กับ 0 ตามลำดับ โมเดลของค่าสังเกต \(y_i\) นี้จึงสามารถกำหนดได้ด้วยการแจกแจงความน่าจะเป็นแบบ Bernoulli กล่าวคือ ความน่าจะเป็นที่จะออกหน้าหัวมีค่าเท่ากับ \(p(y_i=1|\theta)=\theta\) และความน่าจะเป็นที่ออกหน้าก้อยมีค่าเท่ากับ \(p(y_i=0|\theta)=1-\theta\) ซึ่งเขียนเป็นสมการรวมได้ดังนี้

\(p(y_i|\theta)=\theta^{y_i}(1-\theta)^{1-y_i}\) เมื่อ \(i = 1,2,3,...,n\)

จากโมเดลในข้างต้นจะเห็นว่าการเกิดค่าสังเกต \(y_i\) ขึ้นอยู่กับพารามิเตอร์ \(\theta\) ที่มีค่าอยู่บนช่วง \([0,1]\) หากพารามิเตอร์ดังกล่าวมีค่าเท่ากับ \(\theta=0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวและก้อยมีค่าเท่ากัน แต่ถ้าหากพารามิเตอร์ \(\theta>0.5\) นั่นหมายถึงโอกาสในการออกหน้าหัวมีมากกว่า จากความหมายดังกล่าวทำให้สามารถมองได้ว่าพารามิเตอร์ \(\theta\) เป็นพารามิเตอร์ความลำเอียงของเหรียญ หากประมาณพารามิเตอร์ดังกล่าวได้จะสามารถสรุปความลำเอียงของเหรียญได้

เนื่องจากเราไม่ได้เก็บรวบรวมข้อมูลเพียงค่าสังเกตเดียว แต่เราทำการทดลองซ้ำ ๆ จำนวน \(n\) ครั้ง เนื่องจากการทดลองแต่ละครั้งเป็นอิสระซึ่งกันและกัน ทำให้โมเดลของชุดค่าสังเกตหรือเรียกอย่างเป็นทางการว่า ฟังก์ชันภาวะความควรจะเป็นดังนี้

\(p(\bf{y}|\theta)=p(y_1|\theta)p(y_2|\theta)...p(y_n|\theta)\)

ซึ่งมีค่าเท่ากับ

\(p(\bf{y}|\theta)=\theta^{\sum_{i=1}^ny_i}(1-\theta)^{n-\sum_{i=1}^ny_i}\)

ในฟังก์ชันภาวะความควรจะเป็นนี้จะให้ค่าความน่าจะเป็นหรือค่าความหนาแน่นของข้อมูลตัวอย่าง \(\bf{y}\) บนแต่ละค่าที่เป็นไปได้ของพารามิเตอร์ \(\theta\) เนื่องจากข้อมูลตัวอย่างเป็นค่าคงที่ ดังนั้นค่าความหนาแน่นนี้จึงเปลี่ยนแปลงไปตามค่าของพารามิเตอร์ \(\theta\) โดยหากค่าความหนาแน่นนี้มีค่าสูง ณ ค่า \(\theta\) ค่าใด นั่นหมายถึงว่าพารามิเตอร์ค่านั้นทำให้โมเดลของค่าสังเกตกับค่าสังเกตจริงมีความสอดคล้องกันสูง ข้อสังเกตนี้ถูกนำไปใช้เพื่อหาค่าประมาณด้วยวิธีภาวะความควรจะเป็นสูงสุด (maximum likelihood: ML) แต่ในวิธีการแบบเบส์จะนำสารสนเทศส่วนนี้ไปประมวลร่วมกับความเป็นไปได้ของพารามิเตอร์ในการแจกแจงความน่าจะเป็นเบื้องต้นก่อน

กำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์

การกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) สามารถทำได้หลายลักษณะ ทั้งแบบไม่ต่อเนื่อง และแบบต่อเนื่อง ขึ้นอยู่กับขอบเขตของการวิเคราะห์ ข้อมูลในอดีตที่ใช้สนับสนุน และธรรมชาติของพารามิเตอร์นั้น ในตัวอย่างข้างต้นพารามิเตอร์ \(\theta\) ต้องมีค่าอยู่บนช่วง \([0,1]\)

เพื่อให้ตัวอย่างนี้ทำความเข้าใจได้ง่าย จึงเลือกกำหนดการแจกแจงความน่าจะเป็นเบื้องต้นของพารามิเตอร์ \(\theta\) เป็นการแจกแจงแบบไม่ต่อเนื่อง โดยมีค่าที่เป็นไปได้เท่ากับ 0.0, 0.1, 0.2, 0.3, …., 1.0 และมีการแจกแจงของความน่าจะเป็นก่อนหน้าดังรูป

# prior distribution
theta<-seq(0,1,0.1)
p.theta<-c(0.01,0.04,0.075,0.1,0.15,0.25,0.15,0.1,0.075,0.04,0.01)
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab=expression(p(theta)), 
    type="h", lwd=2, col="orange",
    main="Prior Distribution")
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

การประมาณการแจกแจงความน่าจะเป็นภายหลัง

จากการทดลองโยนเหรียญจำนวน 10 ครั้งพบว่า ออกหน้าหัวจำนวน 6 ครั้ง ดังนั้นจะได้ว่า

\(p(\bf{y}|\theta)=\theta^{6}(1-\theta)^4\)

จากขอบเขตของ \(theta\) ที่กำหนดในการแจกแจงความน่าจะเป็นเบื้องต้นจะได้ว่าฟังก์ชันภาวะความควรจะเป็น มีลักษณะดังรูป

# likelihood function
L<-theta^6*(1-theta)^4
plot(theta,L, 
     xlab=expression(theta), 
    ylab="p(y|theta)",
    type="h", lwd=2, col="orange",
    main="Likelihood Function")
points(theta,L, type="p",pch=16, cex=2, col="orange")

จากทฤษฎีบทของเบส์เราสามารถรวมสารสนเทศจากการแจกแจงความน่าจะเป็นเบื้องต้น กับฟังก์ชันภาวะความควรจะเป็นได้ดังนี้

# marginal y
p.y<-sum(p.theta*L)

par(mfrow=c(1,2))
# posterior of theta given y
post.theta<-L*p.theta/p.y
plot(theta,post.theta, 
     xlab=expression(theta), 
    ylab="p(theta|y)",
    type="h", lwd=2, col="blue",
    ylim=c(0,0.5),
    main="Posterior Distribution")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")

# prior and posterior in the same plot
plot(theta,p.theta, 
     xlab=expression(theta), 
    ylab="Probability", 
    type="l", lwd=2, col="orange",
    ylim=c(0,0.5))
points(theta,p.theta, type="p",pch=16, cex=2, col="orange")

points(theta,post.theta, 
    type="l", lwd=2, col="blue")
points(theta,post.theta, type="p",pch=16, cex=2, col="blue")
legend(-0.05,0.48, legend=c("Prior","Posterior"), col=c("orange","blue"), lty=1, bty="n", lwd=2)

การอนุมานเชิงสถิติจากการแจกแจงความน่าจะเป็นภายหลัง

เมื่อได้การแจกแจงความน่าจะเป็นภายหลังที่ต้องการแล้ว ผู้วิเคราะห์สามารถวิเคราะห์การแจกแจงดังกล่าวเพื่ออนุมานเกี่ยวกับพารามิเตอร์ความลำเอียงที่ต้องการศึกษาได้ โดยสามารถทำได้หลายลักษณะ ทั้งการประมาณค่าแบบจุด การประมาณค่าแบบช่วง และการทดสอบสมมุติฐาน รายละเอียดดังนี้

(1) การประมาณค่าแบบจุด

การประมาณค่าแบบจุดสามารถทำได้โดยใช้สถิติพื้นฐานสรุปแนวโน้มสู่ส่วนกลางของค่าพารามิเตอร์ จากการแจกแจงความน่าจะเป็นภายหลัง เช่น ค่าเฉลี่ย มัธยฐาน และฐานนิยม นอกจากนี้ยังสามารถคำนวณค่าส่วนเบี่ยงเบนมาตรฐานเพื่อประเมินระดับความน่าเชื่อถือของค่าประมาณแบบจุดดังกล่าวได้อีกด้วย จากการแจกแจงความน่าจะเป็นภายหลังข้างต้น จะได้ว่า

# mean = expected value of theta
mean.theta<-sum(post.theta*theta)
sd.theta<-sqrt(sum(post.theta*(theta-mean.theta)^2))
mean.theta

## [1] 0.5540441

sd.theta

## [1] 0.1145753

(2) การประมาณค่าแบบช่วง

ช่วงการประมาณที่ได้จากวิธีแบบเบส์จะเรียกว่า ช่วงความน่าเชื่อถือ (credible interval) การประมาณค่าพารามิเตอร์แบบช่วงจากการแจกแจงความน่าจะเป็นภายหลังสามารถทำได้หลายวิธีการ วิธีการง่าย ๆ คือการใช้ค่าเฉลี่ยบวกลบด้วยส่วนเบี่ยงเบนมาตรฐานของพารามิเตอร์ อีกวิธีการหนึ่งที่มีประสิทธิภาพมากกว่าคือการหาช่วงแบบ highest density interval (HDI) กล่าวคือเป็นช่วงการประมาณที่มีความหนาแน่นมากที่สุดที่ทำให้ค่าความน่าจะเป็นของพารามิเตอร์ภายในช่วงดังกล่าวเท่ากับค่าที่กำหนดเช่น ช่วง 95% HDI ของพารามิเตอร์ \(\theta\) ในข้างต้นจะมีค่าประมาณ \([0.4,0.8]\)

ช่วงดังกล่าวมีความแตกต่างจาก 95% confidence interval อย่างไร??

(3) การทดสอบสมมุติฐาน

การทดสอบสมมุติฐานด้วยวิธีการแบบเบส์มีความยืดหยุ่นสูงมาก และสามารถทำได้หลายวิธีการ วิธีการแรกเรียกว่า maximum a posteriori (MAP) test มีรายละเอียดดังนี้

กำหนดให้ \(H_0\) และ \(H_1\) เป็นสมมุติฐานที่ต้องการเปรียบเทียบ ในกรณีนี้จะเห็นว่ามีความไม่แน่นอนว่าสมมุติฐานใดเป็นสมมุติฐานที่ถูกต้องกันแน่ ดังนั้นจึงสามารถกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าของสมมุติฐานทั้งสองได้ ดังนี้

\(p(H_0) = p_0\) และ \(p(H_1)=p_1\) โดยที่ \(p_0+p_1=1\)

และกำหนดให้ \(p(\bf{y}|H_0)\) กับ \(p(\bf{y}|H_1)\) คือฟังก์ชันภาวะความควรจะเป็นบนสมมุติฐานทั้งสอง ซึ่งจากทฤษฎีบทของเบส์จะได้ว่า

\(p(H_0|\bf{y})=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y})}\)

\(p(H_1|\bf{y})=\frac{p(\bf{y}|H_1)p(H_1)}{p(\bf{y})}\)

การตัดสินใจว่าควรเลือกเชื่อสมมุติฐานใด สามารถทำได้โดยดูจากอัตราส่วนที่เรียกว่า posterior odd ซึ่งมีค่าเท่ากับ

\(\frac{p(H_0|\bf{y})}{p(H_1|\bf{y})}=\frac{p(\bf{y}|H_0)p(H_0)}{p(\bf{y}|H_1)p(H_1)}\)

โดยหากอัตราส่วนข้างต้นมีค่ามากกว่า 1 แสดงว่า สมมุติฐาน \(H_0\) มีแนวโน้มน่าเชื่อถือมากกว่า ในทางกลับกันหากอัตราส่วนดังกล่าวมีค่าน้อยกว่า 1 แสดงว่าสมมุติฐาน \(H_1\) มีแนวโน้มน่าเชื่อถือมากกว่า

จากตัวอย่างโยนเหรียญ หากต้องการทดสอบว่าเหรียญมีความเที่ยงตรงหรือไม่ อาจกำหนดสมมุติฐานเป็นดังนี้

\(H_0: \theta = 0.5\) vs \(H_1: \theta>0.5\)

prior.H<-c(0.5,0.5) # prior for H0 and H1 respectively.
# calculate likelihood
L.H0<-0.5^6*(1-0.5)^4
theta.H1<-theta[theta>0.5]
L.H1<-sum(theta.H1^6*(1-theta.H1)^4)

# calculate posterior odd
(prior.H[1]*L.H0)/(prior.H[2]*L.H1)

## [1] 0.3727444

อีกวิธีการหนึ่งที่สามารถทำได้เรียกว่า minimum cost hypothesis test รายละเอียดมีดังนี้

กำหนดให้

• \(C_{10}\) คือ cost ของการเลือก \(H1\) โดยที่ \(H_0\) ถูกต้อง (cost ของการเกิด type I error)
• \(C_{01}\) คือ cost ของการเลือก \(H_0\) โดยที่ \(H_1\) ถูกต้อง (cost ของการเกิด type II error)

เกณฑ์การพิจารณาจะใช้ค่าความเสี่ยงภายหลัง (posterior risk) ดังนี้

\(\frac{p(H_0|\bf{y})C_{10}}{p(H_1|\bf{y})C_{01}}\)

จากตัวอย่างโยนเหรียญสมมุติว่า \(C_{10}=5C_{01}\)

c10<-5
c01<-1

# posterior risk of accepting H0
(prior.H[2]*L.H1)*c01

## [1] 0.001309962

# posterior risk of accepting H1
(prior.H[1]*L.H0)*c10

## [1] 0.002441406

# risk ratio
(prior.H[1]*L.H0)*c10/(prior.H[2]*L.H1)*c01

## [1] 1.863722

การแจกแจงความน่าจะเป็นก่อนหน้า

การแจกแจงความน่าจะเป็นก่อนหน้าใช้สะท้อนความเชื่อเบื้องต้นในพารามิเตอร์ต่าง ๆ ของโมเดล ความเชื่อดังกล่าวมาได้จากทั้งผลการศึกษาในอดีต หลักเหตุผล และประสบการณ์ของนักวิจัย นอกจากนี้ยังควรต้องพิจารณาธรรมชาติหรือค่าที่เป็นไปได้ของพารามิเตอร์ประกอบด้วย

โดยทั่วไปการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าสามารถทำได้ 2 ลักษณะ ลักษณะแรกเรียกว่า การแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ (non-informative prior distribution) และการแจกแจงความน่าจะเป็นก่อนหน้าแบบให้สารสนเทศ (informative prior distribution)

การแจกแจงความน่าจะเป็นก่อนหน้าของพารามิเตอร์อิทธิพลคงที่สามารถเลือกใช้ได้หลายการแจกแจง เช่น การแจกแจงแบบ uniform ที่เป็น noninformative prior การแจกแจงแบบปกติ ที่เป็นไปได้ทั้ง noninformative และ informative prior ขึ้นอยู่กับการกำหนดพารามิเตอร์ค่าเฉลี่ยและส่วนเบี่ยงเบนมาตรฐานของการแจกแจง

จากตัวอย่างข้างต้นกำหนดให้

\(\beta_0 \sim N(0,10000)\)

\(\beta_1 \sim N(0,10000)\)

ซึ่งเป็นการแจกแจงความน่าจะเป็นก่อนหน้าแบบไม่ให้สารสนเทศ

การแจกแจงความน่าจะเป็นก่อนหน้าของส่วนเบี่ยงเบนมาตรฐานหรือความแปรปรวน มีความจำเป็นต้องเลือกการแจกแจงที่ domain ของการแจกแจงมีค่าไม่ติดลบ โดยปกติมักกำหนดให้ความแปรปรวนดังกล่าวมีการแจกแจงแบบแกมมาผกผัน (inverse gamma distribution)

จากตัวอย่างข้างต้น กำหนดให้

\(sigma^2 \sim IG(a,b)\)

โดยที่ \(IG\) คือการแจกแจงแบบ inverse-gamma

Note: ในเชิงทฤษฎีหากกำหนดให้ \(sigma^2 \sim IG(a,b)\) จะได้ว่า

\(p(\sigma^2|a,b) \propto (\sigma^2)^{-a-1}exp(\frac{-b}{\sigma^2})\)

หากกำหนดให้ \(a \rightarrow 0\) และ \(b \rightarrow 0\) แล้วฟังก์ชันความหนาแน่นในข้างต้นจะเขียนใหม่ได้เป็น

\(p(\sigma^2|a,b) \propto \frac{1}{\sigma^2}\) เรียกว่า Jeffery prior ซึ่งเป็น noninformative prior ที่ใช้ได้ตัวหนึ่งของพารามิเตอร์ความแปรปรวน

ดังนั้นในกรณีนี้จึงกำหนดให้

\(sigma^2 \sim IG(0.0001,0.0001)\)

รูปด้านล่างแสดงรายการของการแจกแจงความน่าจะเป็นที่สามารถกำหนดได้ใน BUGS และ JAGS

จากการระบุโมเดลข้างต้นทำให้สามารถเขียน syntax เพื่อระบุโมเดลใน JAGS ได้ดังนี้

model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~gamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}

ขั้นตอนการประมวลผลด้วยโปรแกรม JAGS

การประมวลผลด้วย JAGS มี 4 ขั้นตอนหลักได้แก่

1. เตรียมข้อมูล

2. ระบุโมเดล

ขั้นตอนนี้สามารถทำได้ดังตัวอย่างข้างต้น และเมื่อระบุโมเดลในโปรแกรม R เรียบร้อยแล้ว ผู้วิเคราะห์จำเป็นต้องเขียนโมเดลดังกล่าวในรูปของไฟล์ .txt โดยใช้คำสั่ง writeLines(modelString, con="model.txt") ดังตัวอย่างด้านล่าง

modelString<-"model{
  
  for(i in 1:n)
  {
    # likelihood function
    y[i]~dnorm(mu[i],tau) #where tau is precision = 1/sigma2
    
    # deteministic model
    mu[i]<-b0+b1*x[i] 
  }
  
  # prior distribution
  b0~dnorm(0,10^-4)
  b1~dnorm(0,10^-4)
  
  tau~dgamma(0.0001,0.0001)
  sigma2<-1/tau # deteministic model
}" #Wrap model syntax into String object

writeLines(modelString, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/simReg.txt") # write modelString to model.txt

3. ส่งข้อมูลและโมเดลไปประมวลผลบน JAGS

4. ตรวจสอบการลู่เข้าของ MCMC

5. วิเคราะห์และแปลผล

สมมุติว่าข้อมูลที่ใช้ในการวิเคราะห์เป็นดังรูป โดยเก็บไว้ในตัวแปร dat.reg

head(dat.reg)

##        ATT      ACH
## 1 30.69309 56.96164
## 2 11.95094 28.41848
## 3 32.41264 42.84734
## 4 20.02710 51.75678
## 5 19.52567 32.15973
## 6 37.61670 57.54749

par(mar=c(4,4,1,1))
plot(ATT,ACH, pch=16, xlab="ATT", ylab="ACH")

การประมวลผลโมเดลข้างต้นผ่านโปรแกรม JAGS จะต้องใช้คำสั่ง 2 ตัวได้แก่ฟังก์ชัน jag.model() และ coda.samples() ที่มีอากิวเมนท์สำคัญดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(rjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

model<-jags.model(file = "simReg.txt",
                  data = dataList,
                  inits = initsList,
                  n.chains = 3)
update(model, n.iter=2000) # n.burnin
samples<-coda.samples(model,
                      variable.names=c("b0","b1","sigma2"),
                      n.iter = 10000,
                      thin=1)
head(samples)
plot(samples)
summary(samples)

อีก package หนึ่งที่สามารถเรียก JAGS จาก R ได้คือ package-runjags ซึ่งมีจุดเด่นคือสามารถสั่งประมวลแบบคู่ขนาน (parallel) แยกตาม core ของ CPU ได้ในกรณีที่ผู้วิเคราะห์กำหนดให้สุ่มตัวอย่างจากลูกโซ่มาร์คอฟที่มีจำนวนมากกว่า 1 ลูกโซ่ ตัวอย่างคำสั่งเป็นดังนี้

setwd("/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/")

library(runjags)
dataList<-list(y=dat.reg$ACH, x=dat.reg$ATT, n=length(dat.reg$ACH))
initsList<-list(b0=1, b1=1, tau=1)

runjag.output<-run.jags(method="parallel",
                model="simReg.txt",
                monitor=c("b0","b1","sigma2"),
                data=dataList,
                inits=initsList,
                n.chains=3,
                burnin=2000,
                sample=10000, #n.iter
                thin=1,
                summarise=FALSE,
                plots=FALSE)
  
samples<-as.mcmc.list(runjag.output)
head(samples)
plot(samples)
summary(samples)

การวิเคราะห์ผลลัพธ์ที่ได้จาก MCMC

การวิเคราะห์ผลลัพธ์ที่ได้จาก MCMC มีวัตถุประสงค์ 2 ประการ ประการแรกคือการวิเคราะห์เพื่อตรวจสอบ/วินิจฉัยว่าตัวอย่างที่ได้จาก MCMC ในข้างต้น เป็นตัวอย่างที่สุ่มได้จากการแจกแจงความน่าจะเป็นภายหลังที่เป็นเป้าหมายแล้วหรือไม่ (ลูกโซ่ที่สร้างขึ้นลู่เข้าสู่ posterior distribution เป้าหมายแล้วหรือไม่) การตรวจสอบดังกล่าว

ผลลัพธ์ที่ได้จากอัลกอริทึม MCMC จะเป็นตัวอย่างของพารามิเตอร์ในโมเดลที่มีจำนวนเท่ากับ จำนวนรอบทวนซ้ำ (iteration) ตามที่กำหนด หากจำนวนรอบที่ทวนซ้ำมีมากเพียงพอตัวอย่างดังกล่าวจะเป็นค่าพารามิเตอร์ที่ถูกสุ่มมาจากการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ภายในโมเดล ผู้วิเคราะห์สามารถใข้ข้อมูลตัวอย่างดังกล่าวในการอนุมานเชิงสถิติผ่านการแจกแจงของตัวอย่างพารามิเตอร์ อย่างไรก็ตามก่อนการใช้งานตัวอย่างดังกล่าวผู้วิเคราะห์จำเป็นต้องตรวจสอบคุณสมบัติของตัวอย่างพารามิเตอร์ดังกล่าวก่อน เพื่อรับประกันความถูกต้องของผลการวิเคราะห์ โดยมีประเด็นการตรวจสอบดังนี้

1. อัตสหสัมพันธ์ระหว่างตัวอย่างภายในแต่ละลู่โซ่ (autocorrelation) เนื่องจากตัวอย่างของพารามิเตอร์ถูกจำลองขึ้นโดยสุ่มจากกระบวนการมาร์คอฟ ดังนั้นจึงมีความเป็นไปได้ที่ตัวอย่างในแต่ละรอบจะมีความสัมพันธ์ซึ่งกันและกัน หากความสัมพันธ์ดังกล่าวมีค่าสูงมากเกินไปจะทำให้ผลการวิเคราะห์จากลูกโซ่ดังกล่าวมีความคลาดเคลื่อนได้ การแก้ปัญหาดังกล่าวสามารถทำได้ด้วยการลบตัวอย่าง

และผู้วิเคราะห์สามารถใช้ตัวอย่างของพารามิเตอร์ดังกล่าวแทนค่าประมาณของการแจกแจงความน่าจะเป็นภายหลัง

ตัวอย่างโยนเหรียญ (รอบที่ 2)

• จากตัวอย่างที่ 1 จะได้ว่าโมเดลของค่าสังเกตเป็นฟังก์ชันความน่าจะเป็นแบบ Bernuolli ที่มีพารามิเตอร์เป็น \(\theta\) เขียนแทนด้วย

\(y_i \sim Ber(\theta)\) โดยที่ \(i=1,2,...,n\)

• อย่างไรก็ตามในกรณีนี้ผู้วิเคราะห์ต้องการกำหนดการแจกแจงความน่าจะเป็นก่อนหน้าแบบต่อเนื่องให้กับพารามิเตอร์ความลำเอียง \(\theta\) เนื่องจากพารามิเตอร์ดังกล่าวต้องมีค่าอยู่บนช่วง \([0,1]\) การแจกแจงความน่าจะเป็นที่เหมาะสมจึงเป็นการแจกแจงแบบบีต้า (beta distribution) ที่มีพารามิเตอร์ตำแหน่ง \(a\) และพารามิเตอร์สเกล \(b\) ดังนั้น

\(\theta \sim Beta(a,b)\)

พารามิเตอร์ทั้งสองตัวของการแจกแจงแบบบีต้าทำให้รูปทรงการแจกแจงมีความแตกต่างกัน ซึ่งสามารถใช้สะท้อนความเชื่อเกี่ยวกับความลำเอียงของเหรียญที่แตกต่างกันได้เป็นอย่างดี รูปด้านล่างแสดงตัวอย่างโค้งความหนาแน่นของการแจกแจงแบบบีต้า เมื่อกำหนดพารามิเตอร์ต่าง ๆ

Kruschke (2012)

จากการกำหนดในข้างต้นจะสามารถเขียน syntax ของโมเดลในโปรแรกม JAGS ได้ดังนี้

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~beta(a,b)

# deterministic part
  a<-3
  b<-2
}

การประมวลผลโมเดลข้างต้นด้วยโปรแกรม JAGS บน R จำเป็นต้องมีการ wrap up คำสั่งข้างต้นให้เป็น text file เพื่อส่งออกไปยังโปรแกรม R โดยเขียนคำสั่งดังนี้

coin_model="

model{

# likelihood part    
for (i in 1:n)
{
  y[i]~dbern(theta)
  
}
  
# prior part
  theta~dbeta(a,b)

# deterministic part
  a<-3
  b<-2
}

"
writeLines(coin_model, con="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt")

ขั้นตอนถัดมาคือการประมาณการแจกแจงความน่าจะเป็นภายหลังของพารามิเตอร์ \(\theta\) โดยใช้เทคนิค MCMC ผ่านโปรแกรม JAGS ซึ่งสามารถสั่งงานได้ผ่านฟังก์ชัน jag.model() โดยฟังก์ชันนี้มีอาร์กิวเมนท์ที่สำคัญที่เกี่ยวข้องกับการระบุคุณลักษณะของลูกโซ่มาร์คอฟที่ต้องการสร้างขึ้น ได้แก่

• ข้อมูลค่าสังเกต ในที่นี้สมมุติว่าทำการทดลอง 20 ครั้งได้ผลเป็นดังนี้

y<-rbinom(20,1,0.65)
data.list<-list(y=y, n=20)

• n.chains

library(rjags)
fit<-jags.model(file="/Users/siwachoat/Library/Mobile Documents/com~apple~CloudDocs/markdown/multilevel/modelfiles/coin_model.txt", data=data.list,
                n.chains = 3)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 20
##    Unobserved stochastic nodes: 1
##    Total graph size: 24
## 
## Initializing model

samples<-coda.samples(fit,variable.names=c("theta"),n.iter=5000)

plot(samples)

summary(samples)

## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.7191290      0.0878639      0.0007174      0.0007152 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.5340 0.6610 0.7243 0.7830 0.8750